Apartado 1

El set de datos elegido para este trabajo es el “Goldman Osteometric Data Set”. Este dataset incluye datos relativos a 1538 esqueletos humanos de distintas epocas del Holoceno. Los datos recogidos comprenden variables lineales de 4 huesos largos de las piernas (humero, radio, femur y tibia) y tres variables extra de la pelvis. De manera adicional, se estimo el sexo a partir de la morfologia pelvica y dicha estimacion se incluyo en el Data Set.

El Data Set puede descargarse desde la siguiente pagina: https://web.utk.edu/~auerbach/GOLD.htm#gold

Referencias:

Hemos seleccionado estos datos por el interés científico que consideramos que tienen, además de por su gran valor desde el punto de vista de la Bioestádistica: Esta formado por una muestra poblacional muy extensa, incluye múltiples variables de distinto tipo y puede permitirnos inferir mucha información que antes era desconocida a través de un acercamiento y estudio apropiado.

Apartado 2

El archivo se ha descargado de la dirección citada anteriormente, en formato csv, y se nombra como “data”.

download.file("https://web.utk.edu/~auerbach/Goldman.csv", destfile = "data")

Se carga el archivo en R y se muestran las primeras lineas:

data1 <- read.csv("data")
head(data1)
##   Inst       ID Sex   Age                     NOTE              Location
## 1 AMNH  99.1/69   0 40-50  Tigara - Point Hope, AK Alaska, United States
## 2 AMNH  99.1/70   1   50+  Norton - Point Hope, AK Alaska, United States
## 3 AMNH 99.1/75A   0 30-50 Ipituaq - Point Hope, AK Alaska, United States
## 4 AMNH  99.1/80   0 25-30 Ipituaq - Point Hope, AK Alaska, United States
## 5 AMNH 99.1/83A   1   50+ Ipituaq - Point Hope, AK Alaska, United States
## 6 AMNH 99.1/86B   1 30-40 Ipituaq - Point Hope, AK Alaska, United States
##   Element. LHUM RHUM LRAD RRAD LFEM RFEM LTIB RTIB OSCX Metrics.  LHML LHEB
## 1       NA    0    1    0    1    0    0    1    0    0       NA 308.5   64
## 2       NA    1    1    1    1    0    0    0    0    0       NA    NA   NA
## 3       NA    0    0    1    0    0    0    0    0    1       NA 311.0   59
## 4       NA    0    0    0    0    0    0    0    0    0       NA 289.0   55
## 5       NA    0    0    0    0    0    0    0    0    0       NA 295.0   54
## 6       NA    0    0    0    1    1    0    0    0    0       NA 270.5   55
##    LHHD LHMLD LHAPD RHML RHEB  RHHD RHMLD RHAPD  LRML LRMLD LRAPD  RRML RRMLD
## 1 47.55 24.30 22.00   NA   NA    NA    NA    NA 229.0 14.84 12.02    NA    NA
## 2    NA    NA    NA   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA
## 3 44.44 22.20 22.12  310 59.5 44.11 21.20 22.68    NA    NA    NA 221.5 17.26
## 4 42.94 20.83 20.36  298 58.5 44.41 21.36 22.09 223.5 13.89 11.12 227.0 14.16
## 5 42.51 19.34 19.35  302 55.0 42.06 20.21 19.97 208.0 14.38  9.62 205.0 15.63
## 6 39.74 19.22 19.42  281 55.5 39.84 19.74 19.38 196.5 13.66  8.85    NA    NA
##   RRAPD LFML  LFBL LFEB  LFAB  LFHD LFMLD LFAPD  RFML  RFBL RFEB  RFAB  RFHD
## 1    NA  443 441.5 82.0 70.68 44.92 29.77 31.39 452.0 451.0 85.0 75.85 48.36
## 2    NA  386 383.0 72.0 61.34 44.64 24.97 23.86 383.0 380.5 73.0 62.00 44.09
## 3 11.03  415 410.0 75.0 67.34 43.95 27.64 28.85 415.5 410.0 77.5 67.35 44.17
## 4 11.22  398 392.0 78.5 66.71 43.97 24.82 30.73 400.0 397.0 77.5 66.95 43.77
## 5  9.95  395 393.0 69.0 57.72 40.58 24.69 28.26 396.0 393.0 70.0 57.88 41.13
## 6    NA   NA    NA   NA    NA    NA    NA    NA 375.0 372.5 68.0 58.90 40.80
##   RFMLD RFAPD LTML LTPB LTMLD LTAPD  RTML RTPB RTMLD RTAPD   BIB LIBL RIBL
## 1 28.94 33.01   NA   NA    NA    NA 365.5 77.0 22.41 26.75 274.0  162  161
## 2 24.11 24.41  283   NA 19.12 23.73 281.5   NA 20.04 21.35 266.0   NA   NA
## 3 27.18 28.39  330   73 21.72 26.96 331.5 73.0 21.07 25.36    NA   NA   NA
## 4 25.18 31.50  312   69 21.76 24.56 320.5 71.5 21.98 24.41 240.5   NA  145
## 5 24.79 26.88  306   64 19.44 23.43 305.0   NA 20.30 21.45 285.5  160  155
## 6 23.76 24.40  298   65 17.75 24.64 293.5 62.5 17.05 24.15    NA   NA   NA
##    LAcH  RAcH Derived.  Brachial    Crural  IL.UL.LL IL.LL.UL  CBR.FHD  McH.FHD
## 1 52.56 52.24       NA 0.7423015 0.8167598 0.6611316 1.512558 65.64622 64.52696
## 2 49.21 50.18       NA        NA 0.7340702        NA       NA 65.27654 59.43324
## 3    NA    NA       NA 0.7133655 0.7965081 0.7131367 1.402256 59.28161 58.75034
## 4 50.30 49.80       NA 0.7674617 0.7926065 0.7252709 1.378795 58.81290 58.32493
## 5 45.60 44.83       NA 0.6917923 0.7724399 0.7203994 1.388119 57.61281 51.57435
## 6 45.70 47.04       NA 0.7126020 0.7886667 0.7040626 1.420328 57.49272 51.45120
##   GRINE.FHD  AVG.FHD
## 1  69.27952 66.48423
## 2  64.11982 62.94320
## 3  63.42808 60.48668
## 4  62.99716 60.04500
## 5  56.15914 55.11543
## 6  56.03440 54.99277

Tras ello, se muestra un resumen de las variables y su tipo:

summary(data1)
##      Inst                ID                Sex                Age           
##  Length:1538        Length:1538        Length:1538        Length:1538       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##      NOTE             Location         Element.            LHUM       
##  Length:1538        Length:1538        Mode:logical   Min.   :0.0000  
##  Class :character   Class :character   NA's:1538      1st Qu.:0.0000  
##  Mode  :character   Mode  :character                  Median :0.0000  
##                                                       Mean   :0.1008  
##                                                       3rd Qu.:0.0000  
##                                                       Max.   :1.0000  
##                                                                       
##       RHUM              LRAD             RRAD             LFEM        
##  Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000  
##  Median :0.00000   Median :0.0000   Median :0.0000   Median :0.00000  
##  Mean   :0.08127   Mean   :0.1404   Mean   :0.1268   Mean   :0.07022  
##  3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000  
##  Max.   :1.00000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
##                                                                       
##       RFEM              LTIB              RTIB              OSCX        
##  Min.   :0.00000   Min.   :0.00000   Min.   :0.00000   Min.   :0.00000  
##  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000  
##  Median :0.00000   Median :0.00000   Median :0.00000   Median :0.00000  
##  Mean   :0.06437   Mean   :0.08648   Mean   :0.08648   Mean   :0.02016  
##  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000  
##  Max.   :1.00000   Max.   :1.00000   Max.   :1.00000   Max.   :1.00000  
##                                                                         
##  Metrics.            LHML            LHEB            LHHD           LHMLD      
##  Mode:logical   Min.   :229.5   Min.   :36.51   Min.   :29.58   Min.   :12.01  
##  NA's:1538      1st Qu.:287.4   1st Qu.:53.00   1st Qu.:39.76   1st Qu.:18.23  
##                 Median :303.5   Median :57.50   Median :42.72   Median :19.86  
##                 Mean   :303.8   Mean   :57.44   Mean   :42.74   Mean   :19.88  
##                 3rd Qu.:319.0   3rd Qu.:61.00   3rd Qu.:45.72   3rd Qu.:21.46  
##                 Max.   :376.0   Max.   :75.00   Max.   :56.33   Max.   :27.22  
##                 NA's   :162     NA's   :184     NA's   :170     NA's   :162    
##      LHAPD            RHML            RHEB            RHHD      
##  Min.   :13.00   Min.   :234.0   Min.   :42.00   Min.   :32.24  
##  1st Qu.:18.03   1st Qu.:291.0   1st Qu.:54.00   1st Qu.:39.90  
##  Median :19.64   Median :307.0   Median :58.00   Median :43.02  
##  Mean   :19.72   Mean   :307.6   Mean   :58.19   Mean   :43.00  
##  3rd Qu.:21.30   3rd Qu.:323.5   3rd Qu.:62.00   3rd Qu.:46.00  
##  Max.   :27.29   Max.   :383.0   Max.   :75.00   Max.   :55.67  
##  NA's   :162     NA's   :135     NA's   :154     NA's   :142    
##      RHMLD           RHAPD            LRML           LRMLD      
##  Min.   :12.25   Min.   :13.84   Min.   :179.0   Min.   : 8.74  
##  1st Qu.:18.69   1st Qu.:18.73   1st Qu.:219.0   1st Qu.:12.78  
##  Median :20.30   Median :20.46   Median :233.0   Median :14.05  
##  Mean   :20.36   Mean   :20.46   Mean   :233.1   Mean   :14.10  
##  3rd Qu.:22.04   3rd Qu.:22.19   3rd Qu.:247.0   3rd Qu.:15.32  
##  Max.   :30.44   Max.   :27.19   Max.   :290.5   Max.   :22.15  
##  NA's   :135     NA's   :136     NA's   :217     NA's   :217    
##      LRAPD            RRML         RRMLD           RRAPD            LFML      
##  Min.   : 7.37   Min.   :180   Min.   : 9.16   Min.   : 7.88   Min.   :345.0  
##  1st Qu.:10.23   1st Qu.:221   1st Qu.:13.12   1st Qu.:10.33   1st Qu.:404.5  
##  Median :11.11   Median :235   Median :14.40   Median :11.27   Median :428.0  
##  Mean   :11.19   Mean   :235   Mean   :14.49   Mean   :11.34   Mean   :427.1  
##  3rd Qu.:12.12   3rd Qu.:248   3rd Qu.:15.74   3rd Qu.:12.31   3rd Qu.:448.5  
##  Max.   :16.28   Max.   :291   Max.   :23.22   Max.   :15.68   Max.   :531.0  
##  NA's   :217     NA's   :201   NA's   :198     NA's   :198     NA's   :117    
##       LFBL            LFEB            LFAB            LFHD      
##  Min.   :309.5   Min.   :58.00   Min.   :49.99   Min.   :33.64  
##  1st Qu.:401.0   1st Qu.:72.00   1st Qu.:62.01   1st Qu.:40.61  
##  Median :424.0   Median :76.00   Median :66.61   Median :43.54  
##  Mean   :423.5   Mean   :76.07   Mean   :66.36   Mean   :43.43  
##  3rd Qu.:445.0   3rd Qu.:81.00   3rd Qu.:70.82   3rd Qu.:46.15  
##  Max.   :530.0   Max.   :93.50   Max.   :83.68   Max.   :57.39  
##  NA's   :122     NA's   :158     NA's   :160     NA's   :117    
##      LFMLD           LFAPD            RFML            RFBL      
##  Min.   :17.90   Min.   :18.88   Min.   :341.5   Min.   :279.0  
##  1st Qu.:23.83   1st Qu.:25.12   1st Qu.:403.0   1st Qu.:399.0  
##  Median :25.73   Median :27.16   Median :426.0   Median :422.0  
##  Mean   :25.77   Mean   :27.17   Mean   :425.7   Mean   :421.6  
##  3rd Qu.:27.65   3rd Qu.:29.18   3rd Qu.:447.5   3rd Qu.:443.5  
##  Max.   :34.54   Max.   :37.06   Max.   :532.5   Max.   :531.0  
##  NA's   :115     NA's   :115     NA's   :112     NA's   :115    
##       RFEB            RFAB            RFHD           RFMLD      
##  Min.   :58.00   Min.   :49.83   Min.   :32.97   Min.   :17.97  
##  1st Qu.:72.00   1st Qu.:61.97   1st Qu.:40.51   1st Qu.:23.54  
##  Median :76.00   Median :66.19   Median :43.52   Median :25.37  
##  Mean   :76.26   Mean   :66.28   Mean   :43.50   Mean   :25.44  
##  3rd Qu.:81.00   3rd Qu.:70.61   3rd Qu.:46.40   3rd Qu.:27.24  
##  Max.   :94.00   Max.   :82.85   Max.   :57.86   Max.   :33.84  
##  NA's   :152     NA's   :148     NA's   :103     NA's   :114    
##      RFAPD            LTML            LTPB           LTMLD      
##  Min.   :18.34   Min.   :276.0   Min.   :52.00   Min.   :15.22  
##  1st Qu.:25.10   1st Qu.:333.0   1st Qu.:65.00   1st Qu.:19.78  
##  Median :27.14   Median :353.0   Median :69.50   Median :21.37  
##  Mean   :27.24   Mean   :353.1   Mean   :69.34   Mean   :21.46  
##  3rd Qu.:29.29   3rd Qu.:372.0   3rd Qu.:74.00   3rd Qu.:23.03  
##  Max.   :37.53   Max.   :446.0   Max.   :86.00   Max.   :31.00  
##  NA's   :114     NA's   :135     NA's   :186     NA's   :139    
##      LTAPD            RTML            RTPB           RTMLD      
##  Min.   :19.00   Min.   :237.0   Min.   :50.00   Min.   :15.15  
##  1st Qu.:24.18   1st Qu.:332.0   1st Qu.:65.00   1st Qu.:20.25  
##  Median :26.23   Median :353.0   Median :69.00   Median :21.89  
##  Mean   :26.35   Mean   :352.4   Mean   :69.34   Mean   :21.99  
##  3rd Qu.:28.54   3rd Qu.:372.0   3rd Qu.:74.00   3rd Qu.:23.79  
##  Max.   :37.94   Max.   :444.0   Max.   :85.00   Max.   :29.83  
##  NA's   :140     NA's   :138     NA's   :189     NA's   :143    
##      RTAPD            BIB             LIBL          RIBL            LAcH      
##  Min.   :18.59   Min.   :184.0   Min.   :105   Min.   :106.0   Min.   :38.64  
##  1st Qu.:23.38   1st Qu.:251.0   1st Qu.:145   1st Qu.:145.0   1st Qu.:45.99  
##  Median :25.21   Median :263.0   Median :151   Median :151.0   Median :49.05  
##  Mean   :25.37   Mean   :262.4   Mean   :151   Mean   :151.1   Mean   :48.87  
##  3rd Qu.:27.32   3rd Qu.:274.0   3rd Qu.:158   3rd Qu.:158.0   3rd Qu.:51.72  
##  Max.   :35.20   Max.   :324.0   Max.   :181   Max.   :189.0   Max.   :62.94  
##  NA's   :144     NA's   :69      NA's   :359   NA's   :361     NA's   :167    
##       RAcH       Derived.          Brachial          Crural      
##  Min.   :37.89   Mode:logical   Min.   :0.6779   Min.   :0.6920  
##  1st Qu.:46.05   NA's:1538      1st Qu.:0.7446   1st Qu.:0.8092  
##  Median :49.03                  Median :0.7653   Median :0.8284  
##  Mean   :48.92                  Mean   :0.7653   Mean   :0.8277  
##  3rd Qu.:51.80                  3rd Qu.:0.7857   3rd Qu.:0.8459  
##  Max.   :62.99                  Max.   :1.0769   Max.   :0.9654  
##  NA's   :163                    NA's   :75       NA's   :48      
##     IL.UL.LL         IL.LL.UL        CBR.FHD         McH.FHD     
##  Min.   :0.5993   Min.   :1.299   Min.   :36.09   Min.   :34.69  
##  1st Qu.:0.6803   1st Qu.:1.421   1st Qu.:54.61   1st Qu.:50.99  
##  Median :0.6916   Median :1.446   Median :60.05   Median :57.47  
##  Mean   :0.6922   Mean   :1.446   Mean   :60.15   Mean   :57.41  
##  3rd Qu.:0.7036   3rd Qu.:1.470   3rd Qu.:65.20   3rd Qu.:63.74  
##  Max.   :0.7696   Max.   :1.669   Max.   :92.75   Max.   :89.12  
##  NA's   :120      NA's   :120     NA's   :19      NA's   :19     
##    GRINE.FHD        AVG.FHD     
##  Min.   :39.06   Min.   :38.30  
##  1st Qu.:55.57   1st Qu.:53.74  
##  Median :62.14   Median :59.83  
##  Mean   :62.07   Mean   :59.88  
##  3rd Qu.:68.48   3rd Qu.:65.68  
##  Max.   :94.19   Max.   :92.02  
##  NA's   :19      NA's   :19
str(data1)
## 'data.frame':    1538 obs. of  69 variables:
##  $ Inst     : chr  "AMNH" "AMNH" "AMNH" "AMNH" ...
##  $ ID       : chr  "99.1/69" "99.1/70" "99.1/75A" "99.1/80" ...
##  $ Sex      : chr  "0" "1" "0" "0" ...
##  $ Age      : chr  "40-50" "50+" "30-50" "25-30" ...
##  $ NOTE     : chr  "Tigara - Point Hope, AK" "Norton - Point Hope, AK" "Ipituaq - Point Hope, AK" "Ipituaq - Point Hope, AK" ...
##  $ Location : chr  "Alaska, United States" "Alaska, United States" "Alaska, United States" "Alaska, United States" ...
##  $ Element. : logi  NA NA NA NA NA NA ...
##  $ LHUM     : int  0 1 0 0 0 0 0 0 0 0 ...
##  $ RHUM     : int  1 1 0 0 0 0 0 0 0 0 ...
##  $ LRAD     : int  0 1 1 0 0 0 0 0 0 0 ...
##  $ RRAD     : int  1 1 0 0 0 1 0 0 0 0 ...
##  $ LFEM     : int  0 0 0 0 0 1 0 0 0 0 ...
##  $ RFEM     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ LTIB     : int  1 0 0 0 0 0 0 0 0 0 ...
##  $ RTIB     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ OSCX     : int  0 0 1 0 0 0 1 0 0 0 ...
##  $ Metrics. : logi  NA NA NA NA NA NA ...
##  $ LHML     : num  308 NA 311 289 295 ...
##  $ LHEB     : num  64 NA 59 55 54 55 63 NA 63.5 57 ...
##  $ LHHD     : num  47.5 NA 44.4 42.9 42.5 ...
##  $ LHMLD    : num  24.3 NA 22.2 20.8 19.3 ...
##  $ LHAPD    : num  22 NA 22.1 20.4 19.4 ...
##  $ RHML     : num  NA NA 310 298 302 ...
##  $ RHEB     : num  NA NA 59.5 58.5 55 55.5 60 57 62 58 ...
##  $ RHHD     : num  NA NA 44.1 44.4 42.1 ...
##  $ RHMLD    : num  NA NA 21.2 21.4 20.2 ...
##  $ RHAPD    : num  NA NA 22.7 22.1 20 ...
##  $ LRML     : num  229 NA NA 224 208 ...
##  $ LRMLD    : num  14.8 NA NA 13.9 14.4 ...
##  $ LRAPD    : num  12.02 NA NA 11.12 9.62 ...
##  $ RRML     : num  NA NA 222 227 205 ...
##  $ RRMLD    : num  NA NA 17.3 14.2 15.6 ...
##  $ RRAPD    : num  NA NA 11.03 11.22 9.95 ...
##  $ LFML     : num  443 386 415 398 395 ...
##  $ LFBL     : num  442 383 410 392 393 ...
##  $ LFEB     : num  82 72 75 78.5 69 NA 79 73 79 73 ...
##  $ LFAB     : num  70.7 61.3 67.3 66.7 57.7 ...
##  $ LFHD     : num  44.9 44.6 44 44 40.6 ...
##  $ LFMLD    : num  29.8 25 27.6 24.8 24.7 ...
##  $ LFAPD    : num  31.4 23.9 28.9 30.7 28.3 ...
##  $ RFML     : num  452 383 416 400 396 ...
##  $ RFBL     : num  451 380 410 397 393 ...
##  $ RFEB     : num  85 73 77.5 77.5 70 68 80 74 81 71 ...
##  $ RFAB     : num  75.8 62 67.3 67 57.9 ...
##  $ RFHD     : num  48.4 44.1 44.2 43.8 41.1 ...
##  $ RFMLD    : num  28.9 24.1 27.2 25.2 24.8 ...
##  $ RFAPD    : num  33 24.4 28.4 31.5 26.9 ...
##  $ LTML     : num  NA 283 330 312 306 298 351 361 338 311 ...
##  $ LTPB     : num  NA NA 73 69 64 65 73 69 74.5 68 ...
##  $ LTMLD    : num  NA 19.1 21.7 21.8 19.4 ...
##  $ LTAPD    : num  NA 23.7 27 24.6 23.4 ...
##  $ RTML     : num  366 282 332 320 305 ...
##  $ RTPB     : num  77 NA 73 71.5 NA 62.5 73 70 71 70 ...
##  $ RTMLD    : num  22.4 20 21.1 22 20.3 ...
##  $ RTAPD    : num  26.8 21.4 25.4 24.4 21.4 ...
##  $ BIB      : num  274 266 NA 240 286 ...
##  $ LIBL     : num  162 NA NA NA 160 NA NA 158 160 160 ...
##  $ RIBL     : num  161 NA NA 145 155 NA NA 156 161 160 ...
##  $ LAcH     : num  52.6 49.2 NA 50.3 45.6 ...
##  $ RAcH     : num  52.2 50.2 NA 49.8 44.8 ...
##  $ Derived. : logi  NA NA NA NA NA NA ...
##  $ Brachial : num  0.742 NA 0.713 0.767 0.692 ...
##  $ Crural   : num  0.817 0.734 0.797 0.793 0.772 ...
##  $ IL.UL.LL : num  0.661 NA 0.713 0.725 0.72 ...
##  $ IL.LL.UL : num  1.51 NA 1.4 1.38 1.39 ...
##  $ CBR.FHD  : num  65.6 65.3 59.3 58.8 57.6 ...
##  $ McH.FHD  : num  64.5 59.4 58.8 58.3 51.6 ...
##  $ GRINE.FHD: num  69.3 64.1 63.4 63 56.2 ...
##  $ AVG.FHD  : num  66.5 62.9 60.5 60 55.1 ...

Como puede verse, el dataset tiene una gran cantidad de variables de distintos tipos (1538 observaciones de 69 variables), incluyendo numéricas (las relativas a las mediciones en los huesos), cualitativas (las relativas a la presencia o ausencia de uno de los huesos) y de caracteres (información adicional como el sexo, el código de identificación o el lugar del yacimiento).

La descripcion de las medidas se encuentra en el PDF encontrado en el mismo sitio que los datos: https://web.utk.edu/~auerbach/GoldMeasures.pdf

Podemos ver que las columnas Elements, Metrics y Derived estan compuestas unicamente de NA y sirven unicamente a separar las variables por tipo establecido por los autores (Elements, Metrics y Derived) dentro del dataset.

Otras cosas importantes de notar:

Apartado 3

Preguntas posibles:

  1. ¿Que huesos se preservan mas de la derecha si se preserva la izquierda? (Elements)
  2. ¿Los museos estan especializados en origen de huesos?
  3. (Metrics) la variable BIB (Bi‐iliac Breadth) no tiene derecha ni izquierda. Hay alguna variable con cuya ella esta correlada (que no sea genero)
  4. Variable Sex : podemos hacer una predicion de genero a partir de este dataset para las observaciones donde Sex esta marcado como 0? y 1? ? Los resultados coinciden con el metodo que usaron los autores?
  5. ¿Tienen de media los húmeros derechos e izquierdos (LHML y RHML) la misma longitud máxima o hay diferencias?
  6. ¿Que porcentaje de la población total corresponde a los individuos de +50 años? ¿En qué llacimiento hay mas individuos de este rango de población?

Respuetas:

a) ¿Que huesos se preservan mas de la derecha si se preserva la izquierda? (Elements)

Tenemos que mirar la proporcion de preservacion de las variables de tipo Elements (LHUM, RHUM, LRAD, RRAD, LFEM, RFEM, LTIB, RTIB). La poblacion completa para medir el porcentaje de preservacion de los huesos derechas seria para cada par (LHUM y RHUM, LRAD y RRAD, LFEM y RFEM, LTIB y RTIB) la poblacion donde la preservacion de huesos izquierdos seria completa (== 0)

#Se prepara un nuevo data set, donde esta preservado el LHUM
LHUM_preservado <- data1[data1$LHUM == '0', ]
LHUM_preservado <- LHUM_preservado[!is.na(LHUM_preservado$LHUM),]

# Y aqui preparamos el subdataset de LHUM que solo tiene los RHUM preservados
RHUM_preservado <- LHUM_preservado[LHUM_preservado$RHUM == '0', ]
RHUM_preservado <- RHUM_preservado[!is.na(RHUM_preservado$RHUM),]

# proporcion de RHUM preservados cuando LHUM esta preservado
print("humerus: ")
## [1] "humerus: "
nrow(RHUM_preservado)/nrow(LHUM_preservado)
## [1] 0.934201
#Se prepara un nuevo data set, donde esta preservado el LRAD
LRAD_preservado <- data1[data1$LRAD == '0', ]
LRAD_preservado <- LRAD_preservado[!is.na(LRAD_preservado$LRAD),]

# Y aqui preparamos el subdataset de LRAD que solo tiene los RRAD preservados
RRAD_preservado <- LRAD_preservado[LRAD_preservado$RRAD == '0', ]
RRAD_preservado <- RRAD_preservado[!is.na(RRAD_preservado$RRAD),]

# proporcion de RRAD preservados cuando LRAD esta preservado
print("radius: ")
## [1] "radius: "
nrow(RRAD_preservado)/nrow(LRAD_preservado)
## [1] 0.8933434
#Se prepara un nuevo data set, donde esta preservado el LFEM
LFEM_preservado <- data1[data1$LFEM == '0', ]
LFEM_preservado <- LFEM_preservado[!is.na(LFEM_preservado$LFEM),]

# Y aqui preparamos el subdataset de LFEM que solo tiene los RFEM preservados
RFEM_preservado <- LFEM_preservado[LFEM_preservado$RFEM == '0', ]
RFEM_preservado <- RFEM_preservado[!is.na(RFEM_preservado$RFEM),]

# proporcion de RFEM preservados cuando LFEM esta preservado
print("femur: ")
## [1] "femur: "
nrow(RFEM_preservado)/nrow(LFEM_preservado)
## [1] 0.9426573
#Se prepara un nuevo data set, donde esta preservado el LTIB
LTIB_preservado <- data1[data1$LTIB == '0', ]
LTIB_preservado <- LTIB_preservado[!is.na(LTIB_preservado$LTIB),]

# Y aqui preparamos el subdataset de LTIB que solo tiene los RTIB preservados
RTIB_preservado <- LTIB_preservado[LTIB_preservado$RTIB == '0', ]
RTIB_preservado <- RTIB_preservado[!is.na(RTIB_preservado$RTIB),]

# proporcion de RTIB preservados cuando LTIB esta preservado
print("tibia: ")
## [1] "tibia: "
nrow(RTIB_preservado)/nrow(LTIB_preservado)
## [1] 0.9338078

Vemos que la proporcion mas alta la tienen LFEM y RFEM que son los femures. No deberia ser una sorpresa ya que es el hueso mas largo y fuerte del cuerpo humano asi que lo logico seria que es el que se preserve mas, también proporcialmente entre la derecha y la izquierda.

b) ¿Los museos estan especializados en origen de huesos?

Podemos hacer unos barplots por museo

par(mar=c(5,10,4,4))
barplot(table(data1$Location,data1$Inst),horiz=TRUE,las=2, cex.names=.6,col = rainbow(length(unique(data1$Location))))

El barplot sale muy colorado. Podemos ver que WOAC, SfAP, KU, KSU y CMNH estan muy especializados cada uno en una Location (la barra es mayoritarmente de un unico color). Cual? Para averiguarlo, podemos hacer un subset del dataset componiendolo unicamente de observaciones de estos museos.

inst_woac <- data1[data1$Inst == 'WOAC', ]
inst_sfap <- data1[data1$Inst == 'SfAP', ]
inst_ku <- data1[data1$Inst == 'KU', ]
inst_ksu <- data1[data1$Inst == 'KSU', ]
inst_cmnh <- data1[data1$Inst == 'CMNH', ]

# mostramos la reparticion de los origenes por museo
table(inst_woac$Location,inst_woac$Inst)
##                          
##                           WOAC
##   Kentucky, United States   61
table(inst_sfap$Location,inst_sfap$Inst)
##           
##            SfAP
##   Germany   108
##   Tanzania    1
table(inst_ku$Location,inst_ku$Inst)
##        
##         KU
##   Japan 66
table(inst_ksu$Location,inst_ksu$Inst)
##                      
##                       KSU
##   Ohio, United States  21
table(inst_cmnh$Location,inst_cmnh$Inst)
##                          
##                           CMNH
##   Colorado, United States    3
##   Hawaii, United States      2
##   Ohio, United States      170

Vemos que los museos “especializados” en un origen de los huesos tienden a tener colleciones regionales de la region o pais donde esta el museo. Asi mismo, WOAC – Webb Osteology and Archaeology Collection (que se encuentra en el estado de Kentucky, EEUU) tiene unicamente los huesos con origen en Kentucky. SfAP – Staatssammlung fur Anthropologie und Palaeoanatomie (que es una collecion alemana) los huesos con origen en Alemania (a l’exepcion de un ejemplar de Tanzania), KU – Kyoto University (Kyodai) con origen en Japon, KSU – Kent State University y CMNH – Cleveland Museum of Natural History los huesos con origen de Ohio (que es el estado donde estan).

c) (Metrics) la variable BIB (Bi‐iliac Breadth) no tiene derecha ni izquierda. Hay alguna variable con cuya ella esta correlada (que no sea genero)

#Se importa zoo para realizar la sustitución rápidamente
#install.packages("zoo")

# ponemos las metricas en un dataset a parte
metrics <- data1[,c(18:60)]

library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# reemplazamos los NA en cada variable con la media de la variable
metrics2 <- na.aggregate(metrics)

# miramos la correlacion de cada variable con BIB
cor(metrics2[-39], metrics2$BIB)
##            [,1]
## LHML  0.4711816
## LHEB  0.4696301
## LHHD  0.4859928
## LHMLD 0.4124235
## LHAPD 0.3888506
## RHML  0.4915438
## RHEB  0.4690396
## RHHD  0.5104892
## RHMLD 0.4360081
## RHAPD 0.4171695
## LRML  0.3725394
## LRMLD 0.3487594
## LRAPD 0.3437768
## RRML  0.3657437
## RRMLD 0.3618544
## RRAPD 0.3654109
## LFML  0.4557349
## LFBL  0.4552652
## LFEB  0.5052740
## LFAB  0.4620532
## LFHD  0.5474541
## LFMLD 0.5451732
## LFAPD 0.4323758
## RFML  0.4606696
## RFBL  0.4523228
## RFEB  0.4911563
## RFAB  0.4374224
## RFHD  0.5463900
## RFMLD 0.5285022
## RFAPD 0.4349595
## LTML  0.3842317
## LTPB  0.4499900
## LTMLD 0.3859132
## LTAPD 0.3545380
## RTML  0.3783401
## RTPB  0.4636810
## RTMLD 0.3782224
## RTAPD 0.3675581
## LIBL  0.6501100
## RIBL  0.6679058
## LAcH  0.5372142
## RAcH  0.5320178

La correlacion mas grande de BIB es con RIBL (Right Maximum Iliac Blade Length) (0.6679058) y LIBL (Left Maximum Iliac Blade Length) (0.6501100).

d) Variable Sex : podemos hacer una predicion de genero a partir de este dataset para las observaciones donde Sex esta marcado como 0? y 1? ? Los resultados coinciden con el metodo que usaron los autores?

Podemos usar un algoritmo de clasificacion, por ejempo Support Vector Machine (SVM)

Para aplicar ese algoritmo, hemos seguido en parte las instrucciones encontradas aqui: https://www.geeksforgeeks.org/classifying-data-using-support-vector-machinessvms-in-r/

# separamos el dataset data2 en grupos test y train. El test va a incluir todas las observaciones donde el genero no esta seguro (marcado 0? o 1?) y el train todas las demas observaciones.

data2_test <- data1[data1$Sex == '0?'| data1$Sex == '1?', ]
data2_train <- data1[data1$Sex == '0'| data1$Sex == '1', ]

# en el dataset test reemplazamos los valores 0? y 1? con un el equivalente 0 y 1 para compararlos después de entrenar el algoritmo
data2_test$Sex<-replace(data2_test$Sex, data2_test$Sex=="0?","0") 
data2_test$Sex<-replace(data2_test$Sex, data2_test$Sex=="1?","1") 

# para el dataset final de train y test, necesitamos unicamente las columnas de metricas
metricsgen_train <- data2_train[,c(18:60)]
metricsgen_test <- data2_test[,c(18:60)]

# para la variable que va a ser predecida, tomamos a parte la columna Sex en ambos dataset train y test
sex_train <- data2_train[,c(3)]
sex_test <- data2_test[,c(3)]

# pasamos a tipo factor la variable predecida en ambos train y test
sex_train<-factor(sex_train, levels = c(0, 1))
sex_test<-factor(sex_test, levels = c(0, 1))

# reemplazamos con la media de la variable las observaciones que tienen NA en ambos train y test
library(zoo)
metricsgen_train <- na.aggregate(metricsgen_train)
metricsgen_test <- na.aggregate(metricsgen_test)

head(metricsgen_test)
##      LHML LHEB  LHHD LHMLD LHAPD  RHML RHEB  RHHD RHMLD RHAPD     LRML   LRMLD
## 166 294.0   53 39.17 17.85 18.89 296.0   55 39.48 17.98 19.05 227.5000 13.7800
## 169 290.0   53 37.43 17.85 20.62 294.0   55 36.55 17.65 20.51 226.5000 14.7000
## 171 291.0   56 40.88 20.05 18.71 291.5   59 41.01 21.55 19.61 233.9375 13.8475
## 509 324.5   57 42.05 20.81 18.95 329.0   57 42.15 20.16 18.69 259.0000 15.2400
## 744 293.0   58 36.09 21.05 21.92 300.0   55 40.77 21.34 22.15 233.9375 13.8475
## 994 287.0   50 36.55 16.85 16.64 293.0   52 38.25 18.35 17.25 220.5000 12.6500
##       LRAPD     RRML   RRMLD    RRAPD  LFML     LFBL     LFEB     LFAB     LFHD
## 166 10.4600 236.6875 14.3475 10.83625 416.5 413.0000 72.50000 61.91000 40.75000
## 169 11.1600 230.0000 13.8500 11.38000 411.0 408.5000 69.50000 60.72000 38.54000
## 171 10.7125 230.5000 14.9900 11.04000 403.0 402.0000 71.00000 61.00000 40.46000
## 509 11.6600 263.0000 15.0100 11.31000 460.5 459.0000 80.50000 67.11000 43.27000
## 744 10.7125 236.6875 14.3475 10.83625 424.5 420.6667 73.22222 63.47889 41.23111
## 994 10.2300 228.0000 14.6800 10.74000 411.0 404.5000 72.00000 62.34000 40.44000
##        LFMLD    LFAPD     RFML     RFBL     RFEB     RFAB    RFHD    RFMLD
## 166 24.33000 27.17000 415.9375 412.1875 73.05556 63.83111 40.4575 24.62333
## 169 23.36000 27.77000 410.0000 407.0000 70.00000 60.17000 38.9100 22.57000
## 171 25.56000 26.64000 404.5000 402.0000 73.00000 61.83000 40.7000 24.65000
## 509 25.38000 27.53000 415.9375 412.1875 79.50000 68.73000 40.4575 27.05000
## 744 24.44111 26.79444 402.0000 398.0000 69.00000 60.93000 36.6900 25.98000
## 994 26.94000 29.01000 403.0000 395.0000 70.00000 60.87000 39.7600 25.56000
##        RFAPD  LTML LTPB LTMLD LTAPD     RTML     RTPB    RTMLD    RTAPD   BIB
## 166 27.08889 350.5 64.0 20.62 25.65 357.1111 68.38889 21.41556 24.91889 265.0
## 169 28.18000 359.0 64.5 20.65 25.56 360.5000 65.00000 21.50000 25.39000 266.0
## 171 27.13000 334.0 65.0 22.58 25.67 337.5000 66.00000 22.39000 24.79000 269.5
## 509 26.69000 389.5 76.0 22.00 29.82 391.5000 76.00000 21.40000 29.70000 273.0
## 744 28.19000 371.0 74.0 21.16 24.00 376.0000 68.00000 20.73000 22.79000 257.0
## 994 27.43000 343.0 67.0 19.49 26.14 341.0000 67.50000 20.30000 24.59000 280.5
##         LIBL    RIBL     LAcH  RAcH
## 166 140.0000 147.375 46.56000 46.63
## 169 146.0000 145.000 44.44000 45.69
## 171 143.0000 146.000 47.70000 47.12
## 509 150.0000 151.000 46.79111 49.98
## 744 153.0000 154.000 46.06000 44.87
## 994 145.6667 147.375 46.63000 46.44
head(metricsgen_train)
##       LHML    LHEB     LHHD    LHMLD    LHAPD     RHML     RHEB    RHHD
## 1 308.5000 64.0000 47.55000 24.30000 22.00000 307.5837 58.21747 43.0266
## 2 303.7822 57.4669 42.76761 19.89165 19.72556 307.5837 58.21747 43.0266
## 3 311.0000 59.0000 44.44000 22.20000 22.12000 310.0000 59.50000 44.1100
## 4 289.0000 55.0000 42.94000 20.83000 20.36000 298.0000 58.50000 44.4100
## 5 295.0000 54.0000 42.51000 19.34000 19.35000 302.0000 55.00000 42.0600
## 6 270.5000 55.0000 39.74000 19.22000 19.42000 281.0000 55.50000 39.8400
##      RHMLD    RHAPD     LRML    LRMLD    LRAPD     RRML    RRMLD    RRAPD
## 1 20.36844 20.47182 229.0000 14.84000 12.02000 234.9537 14.49474 11.33899
## 2 20.36844 20.47182 233.0636 14.09854 11.19121 234.9537 14.49474 11.33899
## 3 21.20000 22.68000 233.0636 14.09854 11.19121 221.5000 17.26000 11.03000
## 4 21.36000 22.09000 223.5000 13.89000 11.12000 227.0000 14.16000 11.22000
## 5 20.21000 19.97000 208.0000 14.38000  9.62000 205.0000 15.63000  9.95000
## 6 19.74000 19.38000 196.5000 13.66000  8.85000 234.9537 14.49474 11.33899
##       LFML     LFBL     LFEB     LFAB     LFHD    LFMLD    LFAPD  RFML  RFBL
## 1 443.0000 441.5000 82.00000 70.68000 44.92000 29.77000 31.39000 452.0 451.0
## 2 386.0000 383.0000 72.00000 61.34000 44.64000 24.97000 23.86000 383.0 380.5
## 3 415.0000 410.0000 75.00000 67.34000 43.95000 27.64000 28.85000 415.5 410.0
## 4 398.0000 392.0000 78.50000 66.71000 43.97000 24.82000 30.73000 400.0 397.0
## 5 395.0000 393.0000 69.00000 57.72000 40.58000 24.69000 28.26000 396.0 393.0
## 6 427.1236 423.4733 76.08643 66.38012 43.44486 25.77588 27.17419 375.0 372.5
##   RFEB  RFAB  RFHD RFMLD RFAPD     LTML    LTPB    LTMLD    LTAPD  RTML
## 1 85.0 75.85 48.36 28.94 33.01 353.0617 69.3506 21.46419 26.35661 365.5
## 2 73.0 62.00 44.09 24.11 24.41 283.0000 69.3506 19.12000 23.73000 281.5
## 3 77.5 67.35 44.17 27.18 28.39 330.0000 73.0000 21.72000 26.96000 331.5
## 4 77.5 66.95 43.77 25.18 31.50 312.0000 69.0000 21.76000 24.56000 320.5
## 5 70.0 57.88 41.13 24.79 26.88 306.0000 64.0000 19.44000 23.43000 305.0
## 6 68.0 58.90 40.80 23.76 24.40 298.0000 65.0000 17.75000 24.64000 293.5
##       RTPB RTMLD RTAPD      BIB     LIBL   RIBL     LAcH     RAcH
## 1 77.00000 22.41 26.75 274.0000 162.0000 161.00 52.56000 52.24000
## 2 69.34963 20.04 21.35 266.0000 151.0209 151.13 49.21000 50.18000
## 3 73.00000 21.07 25.36 262.3938 151.0209 151.13 48.88727 48.93736
## 4 71.50000 21.98 24.41 240.5000 151.0209 145.00 50.30000 49.80000
## 5 69.34963 20.30 21.45 285.5000 160.0000 155.00 45.60000 44.83000
## 6 62.50000 17.05 24.15 262.3938 151.0209 151.13 45.70000 47.04000
sex_test
##  [1] 0 0 1 0 1 1 1 0 0 1
## Levels: 0 1
# Feature Scaling: necesitamos escalar las variables en train y test para que esten entre -1 y 1
metricsgen_train = scale(metricsgen_train)
metricsgen_test = scale(metricsgen_test)


# Aplicamos el algoritmo SVM al dataset de train con la variable de prediccion sex_train
#install.packages('e1071')
library(e1071)
 
classifier = svm(formula = sex_train ~ .,
                 data = metricsgen_train,
                 type = 'C-classification',
                 kernel = 'linear')

# predecimos los resultados de test con el classifier creado
y_pred = predict(classifier, newdata = metricsgen_test)

# creamos la matriz de confusion para ver si ha funcionado
cm = table(sex_test, y_pred)
cm
##         y_pred
## sex_test 0 1
##        0 4 1
##        1 2 3

Leyendo a la matriz de confusion, tenemos a 4 verdaderos positivos, 1 falso positivo, 2 falsos negativos y 3 verdaderos negativos. La accuracy es entonces de 7/10 o sea que en 7 casos de los 10, los resultados creados con SVM son los mismos que los que los autores han usado para predecir el genero de las observaciones donde el genero faltaba. Ahora el problema es que no sabiendo si las observaciones marcadas con ? en la columna Sex apartienen de verdad al genero asignado, es dificil saber si nuestro algoritmo a sido eficiente.

Apartado 4

El data set que tenemos contiene muchas variables que podemos cruzar en varios modos para hacernos unas preguntas.

- Unificacion de los formatos de la variable $Age:

Los catagorias de edad no tienen la misma composicion (20‐22; 22‐25; 25‐30; 30‐40; 40‐50; 50+) o sea que entre 20 y 30 anos tenemos a 3 categorias diferentes (20‐22; 22‐25; 25‐30) y todas las edades superiores a 50 estan en una categoria (50+). Mirando las observaciones por categoria, vemos que hay aun mas dispersion de edades:

library(ggplot2)    

barplot(sort(table(data1$Age),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)

Si si usa la variable Age, nos conviene unificar los valores de las observaciones, estableciendo unas categorias fijas, por ejempo <25 (dado que hay alcunas observaciones que estan en la categoria empezando con 18), 25-30, 30-40, 40-50, 50+

data2<-data1
data2$Age<-replace(data2$Age, data2$Age<25,"18-25") 
data2$Age<-replace(data2$Age, data2$Age>=50,"50+") 
data2$Age<-replace(data2$Age, data2$Age>=25 & data2$Age<30,"25-30") 
data2$Age<-replace(data2$Age, data2$Age>=30 & data2$Age<40,"30-40") 
data2$Age<-replace(data2$Age, data2$Age>=40 & data2$Age<50,"40-50") 


#miramos la distribucion de las nuevas categorias de edades ahora
sort(table(data2$Age),decreasing = TRUE)
## 
## 30-40 40-50 25-30   50+ 18-25 
##   576   443   250   157   112
barplot(table(data2$Age), las=2, cex.names=.6)

Tenemos ahora a una distribucion similar a la normal con la mayoria de las observaciones aparteniendo a individuos entre 30 y 50 anos.

- Investigacion de la variable Location

Seria interesante también ver la distribucion de los orígenes de las observaciones (para un usuario no familiarizado con los sitios arqueologicos, en vez de usar la variable Nota que tiene el nombre del sitio arqueologico de origen, usaremos la variable Location)

sort(table(data2$Location),decreasing = TRUE)
## 
##                     Ohio, United States                                 Germany 
##                                     191                                     108 
##                   Alaska, United States                                   Japan 
##                                     104                                      94 
##                                   Egypt                                 Austria 
##                                      81                                      78 
##               New Mexico, United States                 England, United Kingdom 
##                                      73                                      72 
##                 Kentucky, United States             South Dakota, United States 
##                                      61                                      61 
## Aleutian Islands, Alaska, United States                                   Sudan 
##                                      56                                      53 
##                                   Italy                 Illinois, United States 
##                                      51                                      42 
##                                 Belgium                     Utah, United States 
##                                      41                                      35 
##               California, United States                      Philippine Islands 
##                                      31                                      31 
##                               Australia                                    Peru 
##                                      29                                      29 
##                  Arizona, United States                               Argentina 
##                                      23                                      21 
##                                   Chile               New Jersey, United States 
##                                      21                                      17 
##                         Andaman Islands                              Madagascar 
##                                      15                                      14 
##               Washington, United States                Scotland, United Kingdom 
##                                      14                                      13 
##                                  France                                 Ecuador 
##                                      12                                      11 
##                          Canary Islands                 Colorado, United States 
##                                      10                                       8 
##                                  Russia        Democratic Republic of the Congo 
##                                       8                                       6 
##                 Arkansas, United States                         Solomon Islands 
##                                       4                                       4 
##                            South Africa                               Greenland 
##                                       4                                       3 
##                                   China                   Hawaii, United States 
##                                       2                                       2 
##                               Indonesia                                Malaysia 
##                                       1                                       1 
##                        Papua New Guinea                                Tanzania 
##                                       1                                       1 
##                                Tasmania 
##                                       1
par(mar=c(4,10,4,4))
barplot(sort(table(data2$Location),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)

length(unique(data2$Location))
## [1] 45

Vemos que la mayoria de las observaciones tienen origen en EEUU, seguidos por Egipto y Italia.

- Investigacion de la variable $Inst

Finalmente, nos gustaria ver un poco la distribucion de las observaciones por museo (variable Inst)

sort(table(data2$Inst),decreasing = TRUE)
## 
##  NMNH  AMNH  CMNH   MdH    NM  SfAP   NHM MNdAE    KU  WOAC    DC  IRSN   KSU 
##   298   211   175   147   141   109   105    98    66    61    54    52    21
par(mar=c(4,10,4,4))
barplot(sort(table(data2$Inst),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)

Los museos donde se encuentran mas huesos observados son NMNH (Smithsonian National Museum of Natural History), AMNH (American Museum of Natural History), CMNH (Cleveland Museum of Natural History), MdH (Musee de l’Homme) y NM (Naturhistorisches Museum)

- Variable $Sex

sort(table(data2$Sex),decreasing = TRUE)
## 
##   0   1  0?  1? 
## 985 543   5   5
par(mar=c(4,10,4,4))
barplot(sort(table(data2$Sex),decreasing = FALSE),horiz=TRUE, las=2, cex.names=.6)

# porcentaje de observaciones donde no estan seguros del genero
10/1538
## [1] 0.006501951

Vemos que tenemos casi el doble de observaciones provenientes de individuos hombres que de las mujeres. El numero de observaciones donde los autores no estan seguros del genero es inferior a 1% (es de 0.65%). Hemos intentado en apartado 3 determinar el genero con un algoritmo de clasificacion SVM pero no llegamos a los mismos resultados que los autores (teniamos 30% de error)

- Frecuencias de las variables de tipo Elements (LHUM, RHUM, LRAD, RRAD, LFEM, RFEM, LTIB, RTIB, OSCX)

table(data2$LHUM)
## 
##    0    1 
## 1383  155
table(data2$RHUM)
## 
##    0    1 
## 1413  125
table(data2$LRAD)
## 
##    0    1 
## 1322  216
table(data2$RRAD)
## 
##    0    1 
## 1343  195
table(data2$LFEM)
## 
##    0    1 
## 1430  108
table(data2$RFEM)
## 
##    0    1 
## 1439   99
table(data2$LTIB)
## 
##    0    1 
## 1405  133
table(data2$RTIB)
## 
##    0    1 
## 1405  133
table(data2$OSCX)
## 
##    0    1 
## 1507   31
par(mfrow=c(4,2))

barplot(table(data2$LHUM),horiz=TRUE,xlab="LHUM")
barplot(table(data2$RHUM),horiz=TRUE,xlab="RHUM")
barplot(table(data2$LRAD),horiz=TRUE,xlab="LRAD")
barplot(table(data2$RRAD),horiz=TRUE,xlab="RRAD")
barplot(table(data2$LFEM),horiz=TRUE,xlab="LFEM")
barplot(table(data2$RFEM),horiz=TRUE,xlab="RFEM")
barplot(table(data2$LTIB),horiz=TRUE,xlab="LTIB")
barplot(table(data2$RTIB),horiz=TRUE,xlab="RTIB")

barplot(table(data2$OSCX),horiz=TRUE,xlab="OSCX")

Podemos ver que en todas las variables predomina la presencia del elemento y no su ausencia.

En cuanto a las variables del tipo Metrics, podemos ver su distribucion con unos Boxplot. Lo intuitivo puede ser mirar al lado las variables los huesos de izquierda y su equivalente de derecha

# ponemos las metricas en un dataset a parte
metrics <- data2[,c(18:60)]
# reordenamos la metrica para tener la izquierda al lado de la derecha para cada variable
metrics <- metrics[, c(1,6,2,7,3,8,4,9,5,10,11,14,12,15,13,16,17,24,18,25,19,26,20,27,21,28,22,29,23,30,31,35,32,36,33,37,34,38,40,41,42,43,39)]


boxplot(metrics,las=2, cex.names=.2)

Podemos ver ya que hay siempre una cierta disproporcion entre las variables de la izquierda y de la derecha, que sea de tamano promedio o de la varianza.

Apartado 5

Este data set permite múltiples cálculos probabilísticos, dada la gran cantidad de variables y distribuciones que contiene.

Si se selecciona un individuo al azar, ¿Cual es la probabilidad de que su achura bi-ilíaca (BIB) sea mayor que 270 mm?

#Primero se prepara un data set con los datos de la variable BIB (quitando los valores NA):
data_BIB_todaslasvariables <- data1[!is.na(data1$BIB),]
data_BIB <- data_BIB_todaslasvariables$BIB
#Se comprueba normalidad
ks.test(data_BIB, 'pnorm')
## Warning in ks.test(data_BIB, "pnorm"): ties should not be present for the
## Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  data_BIB
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
#como hay normalidad, se calcula media y desviación típica de la distribución para calcular las probabilidades.
mean(data_BIB)
## [1] 262.3958
sd(data_BIB)
## [1] 18.27456
pnorm(q=270, mean = mean(data_BIB), sd = sd(data_BIB), lower.tail = FALSE)
## [1] 0.3386662

La probabilidad es del 33,87%.

Si se selecciona una mujer al azar, ¿Cual es la probabilidad de que su achura bi-ilíaca (BIB) sea mayor que 270 mm?

#Se prepara un nuevo data set, esta vez solo con mujeres, de la variable BIB.
BIB_mujeres <- data1[data1$Sex == '1', ]
data_BIB_variablesmujeres <- BIB_mujeres[!is.na(BIB_mujeres$BIB),]
data_BIB_mujeres <- data_BIB_variablesmujeres$BIB
#Se comprueba normalidad
ks.test(data_BIB_mujeres, 'pnorm')
## Warning in ks.test(data_BIB_mujeres, "pnorm"): ties should not be present for
## the Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  data_BIB_mujeres
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
#como hay normalidad, se calcula media y desviación típica de la distribución para calcular las probabilidades.
mean(data_BIB_mujeres)
## [1] 256.4375
sd(data_BIB_mujeres)
## [1] 17.81262
pnorm(q=270, mean = mean(data_BIB_mujeres), sd = sd(data_BIB_mujeres), lower.tail = FALSE)
## [1] 0.2232096

La probabilidad es del 22,32%.

Si se selecciona un hombre al azar, ¿Cual es la probabilidad de que su achura bi-ilíaca (BIB) esté comprendida entre los 260 y los 280 mm?

#Se prepara otro data set, esta vez solo con hombres, de la variable BIB.
BIB_hombres <- data1[data1$Sex == '0', ]
data_BIB_variableshombres <- BIB_hombres[!is.na(BIB_hombres$BIB),]
data_BIB_hombres <- data_BIB_variableshombres$BIB
#Se comprueba normalidad
ks.test(data_BIB_hombres, 'pnorm')
## Warning in ks.test(data_BIB_hombres, "pnorm"): ties should not be present for
## the Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  data_BIB_hombres
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
#como hay normalidad, se usa la media y desviación típica de la distribución para calcular la probabilidad de interés. Se emplea la fórmula siguiente:
1-pnorm(q=260, mean = mean(data_BIB_hombres), sd = sd(data_BIB_hombres), lower.tail = TRUE)-pnorm(q=280, mean = mean(data_BIB_hombres), sd = sd(data_BIB_hombres), lower.tail = FALSE)
## [1] 0.4160365

La probabilidad es del 41,60%.

Apartado 6

Antes de empezar, para poder realizar las regresiones lineales, es necesario modificar el data frame para excluir los valores NA. Para ello se itera en el data set, sustituyendo cada NA por la media de su columna:

#Se importa zoo para realizar la sustitución rápidamente
library(zoo)
metrics2 <- na.aggregate(metrics)
#Dado que el pool de variables es muy grande, el dataframe metrics2 creado anteriormente se divide en dos.
metrics2_1 <- metrics2[, c(1,6,2,7,3,8,4,9,5,10,11,14,12,15,13,16,17,24,18,25,19,26,20,27)]
metrics2_2 <- metrics2[, c(21,28,22,29,23,30,31,35,32,36,33,37,34,38,40,41,42,43,39)]

#Se observa si hay correlación entre variables.
library(corrplot)
## corrplot 0.92 loaded
corrplot(cor(metrics2_1), method="number",order="hclust")

corrplot(cor(metrics2_2), method="number",order="hclust")

Ademés de la evidente correlación entre lados izquierdo y derecho para las mismas medidas, existe en general una alta correlación entre muchos pares. Como es dificil leer el resultado en una tabla con tantas variables, se hace de nuevo una separación en tablas mas pequeñas.

metrics2_1_1 <- metrics2_1[, c(1,2,3,4,5,6,7,8)]
metrics2_1_2 <- metrics2_1[, c(9,10,11,12,13,14,15,16)]
metrics2_1_3 <- metrics2_1[, c(17,18,19,20,21,22,23,24)]
metrics2_2 <- metrics2[, c(21,28,22,29,23,30,31,35,32,36,33,37,34,38,40,41,42,43,39)]
metrics2_2_1 <- metrics2_2[, c(1,2,3,4,5,6,7,8)]
metrics2_2_2 <- metrics2_2[, c(9,10,11,12,13,14,15,16,17,18,19)]
corrplot(cor(metrics2_1_1), method="number",order="hclust")

corrplot(cor(metrics2_1_2), method="number",order="hclust")

corrplot(cor(metrics2_1_3), method="number",order="hclust")

corrplot(cor(metrics2_2_1), method="number",order="hclust")

corrplot(cor(metrics2_2_2), method="number",order="hclust")

La correlacion mas alta (entre variables que no sean gemelas), se da entre los pares RHEB-RHHD, RFML-RFBL, LFML-LFBL y LFAB-RFEB.

Ahora se realizan las pruebas de regresión lineal entre estos pares:

Regresion1<-lm(RHEB~RHHD, data=metrics2)
summary(Regresion1)
## 
## Call:
## lm(formula = RHEB ~ RHHD, data = metrics2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.260  -1.766   0.000   1.694  10.953 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.8584     0.7948   14.92   <2e-16 ***
## RHHD          1.0775     0.0184   58.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.926 on 1536 degrees of freedom
## Multiple R-squared:  0.6906, Adjusted R-squared:  0.6904 
## F-statistic:  3429 on 1 and 1536 DF,  p-value: < 2.2e-16
plot(RHEB~RHHD, data=metrics2)
abline(11.86, 1.08)
abline(Regresion1)

Regresion2<-lm(RFML~RFBL, data=metrics2)
summary(Regresion2)
## 
## Call:
## lm(formula = RFML ~ RFBL, data = metrics2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -94.007  -1.610  -0.142   1.146  97.314 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.388530   1.747054   5.946 3.39e-09 ***
## RFBL         0.984936   0.004133 238.323  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.963 on 1536 degrees of freedom
## Multiple R-squared:  0.9737, Adjusted R-squared:  0.9737 
## F-statistic: 5.68e+04 on 1 and 1536 DF,  p-value: < 2.2e-16
plot(RFML~RFBL, data=metrics2)
abline(10.39, 0.94)
abline(Regresion2)

Regresion3<-lm(LFML~LFBL, data=metrics2)
summary(Regresion3)
## 
## Call:
## lm(formula = LFML ~ LFBL, data = metrics2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -101.925   -1.327    0.000    1.193   79.624 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.203046   1.506410   5.445 6.01e-08 ***
## LFBL        0.989251   0.003548 278.792  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.218 on 1536 degrees of freedom
## Multiple R-squared:  0.9806, Adjusted R-squared:  0.9806 
## F-statistic: 7.773e+04 on 1 and 1536 DF,  p-value: < 2.2e-16
plot(LFML~LFBL, data=metrics2)
abline(8.20, 0.99)
abline(Regresion3)

Regresion4<-lm(LFAB~RFEB, data=metrics2)
summary(Regresion4)
## 
## Call:
## lm(formula = LFAB ~ RFEB, data = metrics2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.026  -1.499   0.000   1.526  16.209 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.16829    0.94165   5.489 4.73e-08 ***
## RFEB         0.80243    0.01231  65.180  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.858 on 1536 degrees of freedom
## Multiple R-squared:  0.7345, Adjusted R-squared:  0.7343 
## F-statistic:  4248 on 1 and 1536 DF,  p-value: < 2.2e-16
plot(LFAB~RFEB, data=metrics2)
abline(5.17, 0.80)
abline(Regresion4)

Para la regresión múltiple, hemos decidido tratar de predecir la anchura epicondilar del húmero derecho (RHEB) a partir de las otras 7 variables que también han sido usadas en la regresión simple. Se observa el modelo y también se usa el método stepwise para comprobar cuales de las variables son las mejores predictoras.

mrm_metrics2 <- lm(metrics2$RHEB~metrics2$RHHD+metrics2$RFML+metrics2$RFBL+metrics2$LFML+metrics2$LFBL+metrics2$LFAB+metrics2$RFEB)
summary(mrm_metrics2)
## 
## Call:
## lm(formula = metrics2$RHEB ~ metrics2$RHHD + metrics2$RFML + 
##     metrics2$RFBL + metrics2$LFML + metrics2$LFBL + metrics2$LFAB + 
##     metrics2$RFEB)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.7763  -1.7429  -0.0294   1.7024   9.6017 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    3.024e+00  1.055e+00   2.868  0.00419 ** 
## metrics2$RHHD  7.098e-01  3.199e-02  22.186  < 2e-16 ***
## metrics2$RFML  1.490e-03  1.538e-02   0.097  0.92281    
## metrics2$RFBL -1.637e-06  1.421e-02   0.000  0.99991    
## metrics2$LFML -2.025e-02  1.762e-02  -1.150  0.25051    
## metrics2$LFBL  3.018e-02  1.689e-02   1.787  0.07407 .  
## metrics2$LFAB  2.028e-02  2.808e-02   0.722  0.47029    
## metrics2$RFEB  2.431e-01  2.801e-02   8.678  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.754 on 1530 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7258 
## F-statistic: 582.1 on 7 and 1530 DF,  p-value: < 2.2e-16
step(object=mrm_metrics2,direction ="both", trace=1)
## Start:  AIC=3124.18
## metrics2$RHEB ~ metrics2$RHHD + metrics2$RFML + metrics2$RFBL + 
##     metrics2$LFML + metrics2$LFBL + metrics2$LFAB + metrics2$RFEB
## 
##                 Df Sum of Sq   RSS    AIC
## - metrics2$RFBL  1       0.0 11605 3122.2
## - metrics2$RFML  1       0.1 11605 3122.2
## - metrics2$LFAB  1       4.0 11609 3122.7
## - metrics2$LFML  1      10.0 11615 3123.5
## <none>                       11605 3124.2
## - metrics2$LFBL  1      24.2 11629 3125.4
## - metrics2$RFEB  1     571.2 12176 3196.1
## - metrics2$RHHD  1    3733.5 15338 3551.2
## 
## Step:  AIC=3122.18
## metrics2$RHEB ~ metrics2$RHHD + metrics2$RFML + metrics2$LFML + 
##     metrics2$LFBL + metrics2$LFAB + metrics2$RFEB
## 
##                 Df Sum of Sq   RSS    AIC
## - metrics2$RFML  1       0.4 11605 3120.2
## - metrics2$LFAB  1       4.0 11609 3120.7
## - metrics2$LFML  1      10.1 11615 3121.5
## <none>                       11605 3122.2
## - metrics2$LFBL  1      24.4 11629 3123.4
## + metrics2$RFBL  1       0.0 11605 3124.2
## - metrics2$RFEB  1     571.6 12176 3194.1
## - metrics2$RHHD  1    3733.5 15338 3549.2
## 
## Step:  AIC=3120.23
## metrics2$RHEB ~ metrics2$RHHD + metrics2$LFML + metrics2$LFBL + 
##     metrics2$LFAB + metrics2$RFEB
## 
##                 Df Sum of Sq   RSS    AIC
## - metrics2$LFAB  1       3.6 11609 3118.7
## - metrics2$LFML  1       9.9 11615 3119.5
## <none>                       11605 3120.2
## - metrics2$LFBL  1      24.6 11630 3121.5
## + metrics2$RFML  1       0.4 11605 3122.2
## + metrics2$RFBL  1       0.3 11605 3122.2
## - metrics2$RFEB  1     686.5 12292 3206.6
## - metrics2$RHHD  1    3766.0 15371 3550.5
## 
## Step:  AIC=3118.71
## metrics2$RHEB ~ metrics2$RHHD + metrics2$LFML + metrics2$LFBL + 
##     metrics2$RFEB
## 
##                 Df Sum of Sq   RSS    AIC
## - metrics2$LFML  1      10.6 11619 3118.1
## <none>                       11609 3118.7
## + metrics2$LFAB  1       3.6 11605 3120.2
## - metrics2$LFBL  1      26.7 11635 3120.2
## + metrics2$RFBL  1       0.0 11609 3120.7
## + metrics2$RFML  1       0.0 11609 3120.7
## - metrics2$RFEB  1    1074.8 12683 3252.9
## - metrics2$RHHD  1    4083.6 15692 3580.3
## 
## Step:  AIC=3118.11
## metrics2$RHEB ~ metrics2$RHHD + metrics2$LFBL + metrics2$RFEB
## 
##                 Df Sum of Sq   RSS    AIC
## <none>                       11619 3118.1
## + metrics2$LFML  1      10.6 11609 3118.7
## + metrics2$LFAB  1       4.3 11615 3119.5
## + metrics2$RFML  1       1.2 11618 3120.0
## + metrics2$RFBL  1       0.8 11618 3120.0
## - metrics2$LFBL  1     102.7 11722 3129.6
## - metrics2$RFEB  1    1072.5 12692 3251.9
## - metrics2$RHHD  1    4078.2 15697 3578.8
## 
## Call:
## lm(formula = metrics2$RHEB ~ metrics2$RHHD + metrics2$LFBL + 
##     metrics2$RFEB)
## 
## Coefficients:
##   (Intercept)  metrics2$RHHD  metrics2$LFBL  metrics2$RFEB  
##        2.9041         0.7153         0.0120         0.2550

Según el método stepwise, las mejores variables de este conjunto para predecir RHEB son RHHD, LFBL y RFEB.

Apartado 7

Está documentado en Antropología Física desde hace mucho tiempo el hecho de que existen diferencias significativas entre sexos para una gran cantidad de variables óseas lineales. Generalmente, el tamaño de los huesos masculinos es mayor que el de los femeninos. En el data set que hemos empleado hay multitud de ejemplos de este fenómeno, pues suele darse abundantemente en variables de huesos largos.

Por esta razón hemos decidido usar una variable de la cadera (formada por huesos planos en vez de largos), y asi dilucidar si el fenómeno de dimorfismo sexual ocurre también en esta region y en esta población. La variable elegida es la Anchura Bi-ilíaca (BIB).

El primer paso ha sido quitar los casos del data set que tienen los valores “0?” y “1?” para la variable Sex (es decir, los casos cuyo sexo no es seguro).

data_ANOVA<-data1[!(data1$Sex=="0?" | data1$Sex=="1?"),]

Después se ha cambiado la variable Sex de tipo “caracter” a tipo “factor” (necesario para la ANOVA):

data_ANOVA$Sex <- as.factor(data_ANOVA$Sex)

Ahora se puede comenzar con la prueba ANOVA. Hay que recordar que para que la anova pueda implementarse, primero hay que demostrar normalidad entre las dos poblaciones que se comparen (en este caso, hombres y mujeres), asi como igualdad de varianzas en sus distribuciones.

ANOVA

Se crea un data frame nuevo que contenga solo las variables Sex y BIB, a partir del data frame preparado en los pasos anteriores, y que además no contenga casos con NA en dichas variables:

ANOVA_1 <- data.frame(data_ANOVA$Sex, data_ANOVA$BIB)
ANOVA_1 <- na.omit(ANOVA_1)

Se prueba normalidad para ambos sexos en la variable BIB:

library(nortest)
by(data=ANOVA_1,INDICES = ANOVA_1$data_ANOVA.Sex,FUN=function(x){lillie.test(x$data_ANOVA.BIB)})
## ANOVA_1$data_ANOVA.Sex: 0
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$data_ANOVA.BIB
## D = 0.038449, p-value = 0.002299
## 
## ------------------------------------------------------------ 
## ANOVA_1$data_ANOVA.Sex: 1
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x$data_ANOVA.BIB
## D = 0.046259, p-value = 0.00989

La prueba indica que hay normalidad en ambos grupos. Ahora se hace la prueba de igualdad de varianzas:

fligner.test(ANOVA_1$data_ANOVA.BIB~ANOVA_1$data_ANOVA.Sex,ANOVA_1)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  ANOVA_1$data_ANOVA.BIB by ANOVA_1$data_ANOVA.Sex
## Fligner-Killeen:med chi-squared = 0.043118, df = 1, p-value = 0.8355

La prueba indica que hay igualdad de varianzas, por ello, se continua con la ANOVA

anova<-aov(ANOVA_1$data_ANOVA.BIB~ANOVA_1$data_ANOVA.Sex,data=ANOVA_1)
summary(anova)
##                          Df Sum Sq Mean Sq F value Pr(>F)    
## ANOVA_1$data_ANOVA.Sex    1  28664   28664   90.87 <2e-16 ***
## Residuals              1457 459601     315                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El resultado de la ANOVA parece confirmar que existen diferencias significativas en la anchura bi-ilíaca entre hombres y mujeres de esta colección.

Apartado 8

Estos datos han permitido obtener una gran cantidad de resultados y llevar a cabo muchas pruebas de interés gracias a, principalmente, tres factores del data set. El primero es la abundancia en si del registro: la gran cantidad de casos distintos que recoge y la gran variedad de variables tomadas. Esto ha permitido seleccionar las variables mas indicadas para cada prueba que se ha realizado, teniendo siempre disponibles los suficientes individuos para llevarla a cabo. En segundo lugar, las variables que recoge son de distinto tipo, no teniendo únicamente variables continuas como cabe esperar de un data set osteológico, sino tambien de tipo factor (sexo, yacimiento) o numéricas discretas (edad). Esto permite realizar muchas pruebas de alto valor estadístico, como la ANOVA. Por último, dada la naturaleza de los datos que recoge, este data set tiene muchas variables “por parejas”, es decir, que de cada medida ha recogido el valor tanto del lado derecho como del izquierdo. Esto permite realizar estudios de lateralidad, sustituir una variable por su homóloga hermana si no hubiese suficientes casos o se hubiese perdido información, o , si se supiese su lado preferencial (diestro o zurdo) hacer análisis de la influecia de esta variable sobre el desarrollo diferencial de cada lado del esqueleto.

Hemos podido obtener algunas conclusiones intermediarias con el estudio del dataset que hemos hecho. No obstante ello, la cantidad de variables y observaciones presenta muchas mas oportunidades par un estudio mas produndizado que dependeria de la direccion que se quiere tomar (unos huesos especificos? de una origen, genero, tamano especificos? haciendo unos algoritmos predictivos (como en el caso del SVM), usar otro algoritmo (por ejemplo regresion logistica o KNN)? en vez de tener un dataset de train y uno de test, usar también un subset del train para validacion?)